/**
  *  _____  _____  _      __ 
  * |_   _||_   _|| \    /  |
  *   | |    | |  |  \ _/   |
  *   | |    | |  |         |
  *   | |   _| |_ | |\   /| |
  *   \_/   \___/ \_| \_/ |_/
  ******************************************************************************
  * @file    tim_math.c
  * @author  by Tim
  * @brief   fast math algorithm
  * @version V1.0
  ******************************************************************************
  * @attention
  *
  * <h2><center>&copy; Copyright (c) 2021 「Thunder Studio」.
  * All rights reserved.</center></h2>
  *
  *
  *
  *
  ******************************************************************************
  */
  
#include "stdlib.h"
#include "math.h"
#include "stdio.h"



int agaus(double a[],double b[],int n)
{ 
//int n;
//double a[],b[]; //b[]存放解
int *js,l,k,i,j,is,p,q;
double d,t;
	
js=malloc(n*sizeof(int));
l=1;
	
for (i=0;i<n;i++)
	js[i]=i;
	
for (k=0;k<=n-1;k++)
{	 
	d=0.0;
	for (i=k;i<=n-1;i++)
	for (j=k;j<=n-1;j++)
	{ 
		t=fabs(a[i*n+j]);
		if (t>d) 
		{ 
			d=t; 
			p=js[j];
			js[j]=js[k];
			js[k]=p; 
			is=i;
		}
	} //选主元
	
	if (d+1.0==1.0) 
		l=0;
	else
	{ 
		if (js[k]!=k)
			for (i=0;i<=n-1;i++)
			{	
				p=i*n+k; q=i*n+js[k];
				t=a[p]; a[p]=a[q]; a[q]=t;
			} //假若选中的需要列交换，就做列交换
	
		if (is!=k)
		{ 
			for (j=k;j<=n-1;j++)
			{ 
				p=k*n+j; q=is*n+j;
				t=a[p]; a[p]=a[q]; a[q]=t;
			}
			t=b[k]; b[k]=b[is]; b[is]=t;
		} //假若选中的需要行交换，就做行交换
	}
	
	if (l==0)
	{ 
		free(js); 
		printf("fail\n");
		return(0);
	} // 没有唯一解
	
	d=a[k*n+k];
	for (j=k+1;j<=n-1;j++)
	{ 
		p=k*n+j; a[p]=a[p]/d;}
		b[k]=b[k]/d;
		for (i=k+1;i<=n-1;i++)
		{ 
			for (j=k+1;j<=n-1;j++)
			{ 
				p=i*n+j;
				a[p]=a[p]-a[i*n+k]*a[k*n+j];
			}
			b[i]=b[i]-a[i*n+k]*b[k];
		} // 将下述列的元素变为0
	} //全选主元，形成上三角
	
	
	d=a[(n-1)*n+n-1];
	if (fabs(d)+1.0==1.0)
	{ 
		free(js); 
		printf("fail\n");
		return(0);
	}
	
	
	b[n-1]=b[n-1]/d;
	for (i=n-2;i>=0;i--)
	{ 
		t=0.0;
		for (j=i+1;j<=n-1;j++)
			t=t+a[i*n+j]*b[j];
			b[i]=b[i]-t;
	} // 求出解
	
	js[n-1]=n-1;
	for (k=n-1;k>=0;k--)
		if (js[k]!=k)
		{ 
			t=b[k]; 
			b[k]=b[js[k]]; 
			b[js[k]]=t;
		} //还原解
	free(js);
	
	return(1);
}

